{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# relax_box calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The relax_box calculation style refines the lattice parameters of an orthogonal system (crystal structure) by relaxing the box dimensions towards a given pressure.  In refining the lattice parameter values, the box dimensions are allowed to relax, but the relative positions of the atoms within the box are held fixed.\n",
    "\n",
    "This calculations provides a quick tool for obtaining lattice parameters for ideal crystal structures.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Description updated and small changes due to iprPy version.\n",
    "- 2020-05-22: Version 0.10 update - potentials now loaded from database.\n",
    "- 2020-09-22: Setup and parameter definition streamlined.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- With this method there is no guarantee that the resulting parameters are for a stable structure.  Allowing internal relaxations may result in different values for some structures.  Additionally, some transformation paths may be restricted from occurring due to symmetry, i.e. initially cubic structures may remain cubic instead of relaxing to a non-cubic structure.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "The math in this section uses Voigt notation, where indicies i,j correspond to 1=xx, 2=yy, 3=zz, 4=yz, 5=xz, and 6=xy, and x, y and z are orthogonal box vectors.\n",
    "\n",
    "An initial system (and corresponding unit cell system) is supplied with box dimensions, $a_i^0$, close to the equilibrium values.  A LAMMPS simulation is performed that evaluates the system's pressures, $P_{i}$, for the initial system as given, and subjected to twelve different strain states corresponding to one of $\\epsilon_{i}$ being given a value of $\\frac{\\Delta \\epsilon}{2}$, where $\\Delta \\epsilon$ is the strain range parameter. Using the $P_{i}$ values obtained from the strained states, the $C_{ij}$ matrix for the system is estimated as\n",
    "\n",
    "$$ C_{ij} \\approx - \\frac{P_i(\\epsilon_j=\\frac{\\Delta \\epsilon}{2}) - P_i(\\epsilon_j=-\\frac{\\Delta \\epsilon}{2})}{\\Delta \\epsilon}.$$\n",
    "\n",
    "The negative out front comes from the fact that the system-wide stress state is $\\sigma_i = -P_i$. Using $C_{ij}$, an attempt is made to compute the elastic compliance matrix as $S_{ij} = C_{ij}^{-1}$. If successful, new box dimensions are estimated using $S_{ij}$, $a_i^0$, and $P_i$ for the unstrained system\n",
    "\n",
    "$$ a_i = \\frac{a_i^0}{1 - (\\sum_{j=1}^3{S_{ij} P_j})}.$$\n",
    "\n",
    "The system is updated using the new box dimensions. The process is repeated until either $a_i$ converge less than a specified tolerance, $a_i$ diverge from $a_i^0$ greater than some limit, or convergence is not reached after 100 iterations. If the calculation is successful, the final $a_i$ dimensions are reported.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import shutil\n",
    "import datetime\n",
    "from copy import deepcopy\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np  \n",
    "\n",
    "# https://github.com/usnistgov/DataModelDict \n",
    "from DataModelDict import DataModelDict as DM\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'relax_box'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "library = iprPy.Library(load=['lammps_potentials'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use (required).\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is generated using the load parameters and symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.500,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.500,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.500]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.750 |   1.750\n",
      "      2 |       1 |   1.750 |   0.000 |   1.750\n",
      "      3 |       1 |   1.750 |   1.750 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('prototype', 'A1--Cu--fcc', symbols='Ni', a=3.5)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Modify system\n",
    "\n",
    "- __sizemults__ list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "\n",
    "- __system__ is an atomman.System to perform the scan on (required). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of atoms in system = 108\n"
     ]
    }
   ],
   "source": [
    "sizemults = [3, 3, 3]\n",
    "\n",
    "# Generate system by supersizing ucell\n",
    "system = ucell.supersize(*sizemults)\n",
    "print('# of atoms in system =', system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Specify calculation-specific run parameters\n",
    "\n",
    "- __strainrange__ specifies the $\\Delta \\epsilon$ strain range to use in estimating $C_{ij}$.\n",
    "\n",
    "- __pressure_xx__ gives the xx component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __pressure_yy__ gives the yy component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __pressure_zz__ gives the zz component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __convergence_tol__ is the relative tolerance to use in identifying if the lattice constants have converged. \n",
    "\n",
    "- __divergence_scale__ is a factor for identifying if the lattice constants have diverged from the original guess. Divergence is identified if $a > a^0 d$ or $a < a^0 / d$, where d is divergence_scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "strainrange = 1e-7\n",
    "pressure_xx = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_yy = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_zz = uc.set_in_units(0.0, 'GPa')\n",
    "convergence_tol = 1e-11\n",
    "divergence_scale = 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. cij.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('cij.template', 'w') as f:\n",
    "    f.write(\"\"\"# LAMMPS input script that evaluates the Virial pressure of a system both as\n",
    "# it was given, and with small strains applied.\n",
    "# This can be used to estimate the 0 K elastic constants of the system.\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "change_box all triclinic\n",
    "\n",
    "# Specify strain range and number of points to examine in the range\n",
    "variable delta equal <delta>\n",
    "variable steps equal <steps>\n",
    "variable deltax equal ${delta}/(${steps}-1)\n",
    "\n",
    "# Specify variables of the initial configuration's dimensions\n",
    "variable lx0 equal $(lx)\n",
    "variable ly0 equal $(ly)\n",
    "variable lz0 equal $(lz)\n",
    "\n",
    "# Specify the thermo properties to calculate\n",
    "variable peatom equal pe/atoms\n",
    "thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "# Compute properties for the initial configuration\n",
    "run 0\n",
    "\n",
    "# Compute properties for normal x-direction strains\n",
    "variable aratio equal 1-${delta}/2.+(v_a-1)*${deltax}\n",
    "variable xmax equal v_aratio*${lx0}\n",
    "label loopa\n",
    "variable a loop ${steps}\n",
    "change_box all x final 0 ${xmax} remap units box\n",
    "run 0\n",
    "next a\n",
    "jump cij.in loopa\n",
    "change_box all x final 0 ${lx0} remap units box\n",
    "\n",
    "# Compute properties for normal y-direction strains\n",
    "variable bratio equal 1-${delta}/2.+(v_b-1)*${deltax}\n",
    "variable ymax equal v_bratio*${ly0}\n",
    "label loopb\n",
    "variable b loop ${steps}\n",
    "change_box all y final 0 ${ymax} remap units box\n",
    "run 0\n",
    "next b\n",
    "jump cij.in loopb\n",
    "change_box all y final 0 ${ly0} remap units box\n",
    "\n",
    "# Compute properties for normal z-direction strains\n",
    "variable cratio equal 1-${delta}/2.+(v_c-1)*${deltax}\n",
    "variable zmax equal v_cratio*${lz0}\n",
    "label loopc\n",
    "variable c loop ${steps}\n",
    "change_box all z final 0 ${zmax} remap units box\n",
    "run 0\n",
    "next c\n",
    "jump cij.in loopc\n",
    "change_box all z final 0 ${lz0} remap units box\n",
    "\n",
    "# Compute properties for yz shear strains\n",
    "variable eyz equal (-${delta}/2.+(v_d-1)*${deltax})*${lz0}\n",
    "label loopd\n",
    "variable d loop ${steps}\n",
    "change_box all yz final ${eyz} remap units box\n",
    "run 0\n",
    "next d\n",
    "jump cij.in loopd\n",
    "change_box all yz final 0 remap units box\n",
    "\n",
    "# Compute properties for xz shear strains\n",
    "variable exz equal (-${delta}/2.+(v_e-1)*${deltax})*${lz0}\n",
    "label loope\n",
    "variable e loop ${steps}\n",
    "change_box all xz final ${exz} remap units box\n",
    "run 0\n",
    "next e\n",
    "jump cij.in loope\n",
    "change_box all xz final 0 remap units box\n",
    "\n",
    "# Compute properties for xy shear strains\n",
    "variable exy equal (-${delta}/2.+(v_f-1)*${deltax})*${ly0}\n",
    "label loopf\n",
    "variable f loop ${steps}\n",
    "change_box all xy final ${exy} remap units box\n",
    "run 0\n",
    "next f\n",
    "jump cij.in loopf\n",
    "change_box all xy final 0 remap units box\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. relax_box()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relax_box(lammps_command, system, potential,\n",
    "              mpi_command=None, strainrange=1e-6,\n",
    "              p_xx=0.0, p_yy=0.0, p_zz=0.0, p_xy=0.0, p_xz=0.0, p_yz=0.0,\n",
    "              tol=1e-10, diverge_scale=3.):\n",
    "    \"\"\"\n",
    "    Quickly refines static orthorhombic system by evaluating the elastic\n",
    "    constants and the virial pressure.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    strainrange : float, optional\n",
    "        The small strain value to apply when calculating the elastic\n",
    "        constants (default is 1e-6).\n",
    "    p_xx : float, optional\n",
    "        The value to relax the x tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_yy : float, optional\n",
    "        The value to relax the y tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_zz : float, optional\n",
    "        The value to relax the z tensile pressure component to (default is\n",
    "        0.0).\n",
    "    tol : float, optional\n",
    "        The relative tolerance used to determine if the lattice constants have\n",
    "        converged (default is 1e-10).\n",
    "    diverge_scale : float, optional\n",
    "        Factor to identify if the system's dimensions have diverged.  Divergence\n",
    "        is identified if either any current box dimension is greater than the\n",
    "        original dimension multiplied by diverge_scale, or if any current box\n",
    "        dimension is less than the original dimension divided by diverge_scale.\n",
    "        (Default is 3.0).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'a_lat'** (*float*) - The relaxed a lattice constant.\n",
    "        - **'b_lat'** (*float*) - The relaxed b lattice constant.\n",
    "        - **'c_lat'** (*float*) - The relaxed c lattice constant.\n",
    "        - **'alpha_lat'** (*float*) - The alpha lattice angle.\n",
    "        - **'beta_lat'** (*float*) - The beta lattice angle.\n",
    "        - **'gamma_lat'** (*float*) - The gamma lattice angle.\n",
    "        - **'E_coh'** (*float*) - The cohesive energy of the relaxed system.\n",
    "        - **'stress'** (*numpy.array*) - The measured stress state of the\n",
    "          relaxed system.\n",
    "        - **'C_elastic'** (*atomman.ElasticConstants*) - The relaxed system's\n",
    "          elastic constants.\n",
    "        - **'system_relaxed'** (*atomman.System*) - The relaxed system.\n",
    "    \n",
    "    Raises\n",
    "    ------\n",
    "    RuntimeError\n",
    "        If system diverges or no convergence reached after 100 cycles.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Flag for if values have converged\n",
    "    converged = False\n",
    "    \n",
    "    # Define current and old systems\n",
    "    system_current = deepcopy(system)\n",
    "    system_old = None\n",
    "    \n",
    "    system.dump('atom_dump', f='initial.dump')\n",
    "    \n",
    "    for cycle in range(100):\n",
    "        \n",
    "        # Run LAMMPS and evaluate results based on system_old\n",
    "        results = calc_cij(lammps_command, system_current, potential,\n",
    "                           mpi_command=mpi_command,\n",
    "                           p_xx=p_xx, p_yy=p_yy, p_zz=p_zz,\n",
    "                           strainrange=strainrange, cycle=cycle)\n",
    "        system_new = results['system_new']\n",
    "        \n",
    "        # Compare new and current to test for convergence\n",
    "        if np.allclose(system_new.box.vects,\n",
    "                       system_current.box.vects,\n",
    "                       rtol=tol, atol=0):\n",
    "            converged = True\n",
    "            break\n",
    "        \n",
    "        # Compare old and new to test for double-value convergence\n",
    "        elif system_old is not None and np.allclose(system_new.box.vects,\n",
    "                                                    system_old.box.vects,\n",
    "                                                    rtol=tol, atol=0):\n",
    "            # Update current to average of old and new\n",
    "            system_current.box_set(a = (system_new.box.a+system_old.box.a) / 2.,\n",
    "                                   b = (system_new.box.b+system_old.box.b) / 2.,\n",
    "                                   c = (system_new.box.c+system_old.box.c) / 2.,\n",
    "                                   scale=True)\n",
    "            # Calculate Cij for the averaged system\n",
    "            results = calc_cij(lammps_command, system_current, potential,\n",
    "                               mpi_command=mpi_command,\n",
    "                               p_xx=p_xx, p_yy=p_yy, p_zz=p_zz, \n",
    "                               strainrange=strainrange, cycle=cycle+1)\n",
    "            system_new = results['system_new']\n",
    "            converged = True\n",
    "            break\n",
    "        \n",
    "        # Test for divergence\n",
    "        elif system_new.box.a < system.box.a / diverge_scale:\n",
    "            raise RuntimeError('Divergence of box dimensions')\n",
    "        elif system_new.box.a > system.box.a * diverge_scale:\n",
    "            raise RuntimeError('Divergence of box dimensions')\n",
    "        elif system_new.box.b < system.box.b / diverge_scale:\n",
    "            raise RuntimeError('Divergence of box dimensions')\n",
    "        elif system_new.box.b > system.box.b * diverge_scale:\n",
    "            raise RuntimeError('Divergence of box dimensions')\n",
    "        elif system_new.box.c < system.box.c / diverge_scale:\n",
    "            raise RuntimeError('Divergence of box dimensions')\n",
    "        elif system_new.box.c > system.box.c * diverge_scale:\n",
    "            raise RuntimeError('Divergence of box dimensions')\n",
    "        elif results['E_coh'] == 0.0:\n",
    "            raise RuntimeError('Divergence: cohesive energy is 0')\n",
    "        \n",
    "        # If not converged or diverged, current -> old and new -> current\n",
    "        else:\n",
    "            system_old, system_current = system_current, system_new\n",
    "    \n",
    "    # Return values when converged\n",
    "    if converged:\n",
    "        system_new.dump('atom_dump', f='final.dump')\n",
    "        \n",
    "        # Build results_dict\n",
    "        results_dict = {}\n",
    "        results_dict['dumpfile_initial'] = 'initial.dump'\n",
    "        results_dict['symbols_initial'] = system.symbols\n",
    "        results_dict['dumpfile_final'] = 'final.dump'\n",
    "        results_dict['symbols_final'] = system.symbols\n",
    "        \n",
    "        results_dict['lx'] = system_new.box.lx\n",
    "        results_dict['ly'] = system_new.box.ly\n",
    "        results_dict['lz'] = system_new.box.lz\n",
    "        results_dict['xy'] = system_new.box.xy\n",
    "        results_dict['xz'] = system_new.box.xz\n",
    "        results_dict['yz'] = system_new.box.yz\n",
    "        \n",
    "        results_dict['E_coh'] = results['E_coh']\n",
    "        results_dict['measured_pxx'] = results['measured_pxx']\n",
    "        results_dict['measured_pyy'] = results['measured_pyy']\n",
    "        results_dict['measured_pzz'] = results['measured_pzz']\n",
    "        results_dict['measured_pxy'] = results['measured_pxy']\n",
    "        results_dict['measured_pxz'] = results['measured_pxz']\n",
    "        results_dict['measured_pyz'] = results['measured_pyz']\n",
    "        \n",
    "        return results_dict\n",
    "    else:\n",
    "        raise RuntimeError('Failed to converge after 100 cycles')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3. calc_cij()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_cij(lammps_command, system, potential,\n",
    "             mpi_command=None, p_xx=0.0, p_yy=0.0, p_zz=0.0,\n",
    "             strainrange=1e-6, cycle=0):\n",
    "    \"\"\"\n",
    "    Runs cij.in LAMMPS script to evaluate Cij, and E_coh of the current system,\n",
    "    and define a new system with updated box dimensions to test.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    strainrange : float, optional\n",
    "        The small strain value to apply when calculating the elastic\n",
    "        constants (default is 1e-6).\n",
    "    p_xx : float, optional\n",
    "        The value to relax the x tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_yy : float, optional\n",
    "        The value to relax the y tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_zz : float, optional\n",
    "        The value to relax the z tensile pressure component to (default is\n",
    "        0.0).\n",
    "    cycle : int, optional\n",
    "        Indicates the iteration cycle of quick_a_Cij().  This is used to\n",
    "        uniquely save the LAMMPS input and output files.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'E_coh'** (*float*) - The cohesive energy of the supplied system.\n",
    "        - **'stress'** (*numpy.array*) - The measured stress state of the\n",
    "          supplied system.\n",
    "        - **'C_elastic'** (*atomman.ElasticConstants*) - The supplied system's\n",
    "          elastic constants.\n",
    "        - **'system_new'** (*atomman.System*) - System with updated box\n",
    "          dimensions.\n",
    "    \n",
    "    Raises\n",
    "    ------\n",
    "    RuntimeError\n",
    "        If any of the new box dimensions are less than zero.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "\n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data', f='init.dat',\n",
    "                              potential=potential,\n",
    "                              return_pair_info=True)\n",
    "    lammps_variables['atomman_system_pair_info'] = system_info\n",
    "    \n",
    "    lammps_variables['delta'] = strainrange\n",
    "    lammps_variables['steps'] = 2\n",
    "    \n",
    "    # Write lammps input script\n",
    "    template_file = 'cij.template'\n",
    "    lammps_script = 'cij.in'\n",
    "    template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables,\n",
    "                                         '<', '>'))\n",
    "    \n",
    "    # Run lammps\n",
    "    output = lmp.run(lammps_command, lammps_script, mpi_command=mpi_command,\n",
    "                     return_style='model')\n",
    "    shutil.move('log.lammps', 'cij-'+str(cycle)+'-log.lammps')\n",
    "    \n",
    "    # Extract LAMMPS thermo data. Each term ranges i=0-12 where i=0 is undeformed\n",
    "    # The remaining values are for -/+ strain pairs in the six unique directions\n",
    "    lx = uc.set_in_units(np.array(output.finds('Lx')), lammps_units['length'])\n",
    "    ly = uc.set_in_units(np.array(output.finds('Ly')), lammps_units['length'])\n",
    "    lz = uc.set_in_units(np.array(output.finds('Lz')), lammps_units['length'])\n",
    "    xy = uc.set_in_units(np.array(output.finds('Xy')), lammps_units['length'])\n",
    "    xz = uc.set_in_units(np.array(output.finds('Xz')), lammps_units['length'])\n",
    "    yz = uc.set_in_units(np.array(output.finds('Yz')), lammps_units['length'])\n",
    "    \n",
    "    pxx = uc.set_in_units(np.array(output.finds('Pxx')), lammps_units['pressure'])\n",
    "    pyy = uc.set_in_units(np.array(output.finds('Pyy')), lammps_units['pressure'])\n",
    "    pzz = uc.set_in_units(np.array(output.finds('Pzz')), lammps_units['pressure'])\n",
    "    pxy = uc.set_in_units(np.array(output.finds('Pxy')), lammps_units['pressure'])\n",
    "    pxz = uc.set_in_units(np.array(output.finds('Pxz')), lammps_units['pressure'])\n",
    "    pyz = uc.set_in_units(np.array(output.finds('Pyz')), lammps_units['pressure'])\n",
    "    \n",
    "    pe = uc.set_in_units(np.array(output.finds('PotEng')) / system.natoms,\n",
    "                         lammps_units['energy'])\n",
    "    \n",
    "    # Set the six non-zero strain values\n",
    "    strains = np.array([ (lx[2] -  lx[1])  / lx[0],\n",
    "                         (ly[4] -  ly[3])  / ly[0],\n",
    "                         (lz[6] -  lz[5])  / lz[0],\n",
    "                         (yz[8] -  yz[7])  / lz[0],\n",
    "                         (xz[10] - xz[9])  / lz[0],\n",
    "                         (xy[12] - xy[11]) / ly[0] ])\n",
    "    \n",
    "    # Calculate cij using stress changes associated with each non-zero strain\n",
    "    cij = np.empty((6,6))\n",
    "    for i in range(6):\n",
    "        delta_stress = np.array([ pxx[2*i+1]-pxx[2*i+2],\n",
    "                                  pyy[2*i+1]-pyy[2*i+2],\n",
    "                                  pzz[2*i+1]-pzz[2*i+2],\n",
    "                                  pyz[2*i+1]-pyz[2*i+2],\n",
    "                                  pxz[2*i+1]-pxz[2*i+2],\n",
    "                                  pxy[2*i+1]-pxy[2*i+2] ])\n",
    "        \n",
    "        cij[i] = delta_stress / strains[i] \n",
    "    \n",
    "    for i in range(6):\n",
    "        for j in range(i):\n",
    "            cij[i,j] = cij[j,i] = (cij[i,j] + cij[j,i]) / 2\n",
    "    \n",
    "    C = am.ElasticConstants(Cij=cij)\n",
    "    \n",
    "    S = C.Sij\n",
    "    \n",
    "    # Extract the current stress state\n",
    "    stress = -1 * np.array([[pxx[0], pxy[0], pxz[0]],\n",
    "                            [pxy[0], pyy[0], pyz[0]],\n",
    "                            [pxz[0], pyz[0], pzz[0]]])\n",
    "    \n",
    "    s_xx = stress[0,0] + p_xx\n",
    "    s_yy = stress[1,1] + p_yy\n",
    "    s_zz = stress[2,2] + p_zz\n",
    "    \n",
    "    new_a = system.box.a / (S[0,0]*s_xx + S[0,1]*s_yy + S[0,2]*s_zz + 1)\n",
    "    new_b = system.box.b / (S[1,0]*s_xx + S[1,1]*s_yy + S[1,2]*s_zz + 1)\n",
    "    new_c = system.box.c / (S[2,0]*s_xx + S[2,1]*s_yy + S[2,2]*s_zz + 1)\n",
    "    \n",
    "    if new_a <= 0 or new_b <= 0 or new_c <=0:\n",
    "        raise RuntimeError('Divergence of box dimensions to <= 0')\n",
    "    \n",
    "    system_new = deepcopy(system)\n",
    "    system_new.box_set(a=new_a, b=new_b, c=new_c, scale=True)\n",
    "    \n",
    "    results_dict = {}\n",
    "    results_dict['E_coh'] = pe[0]\n",
    "    results_dict['system_new'] = system_new\n",
    "    results_dict['measured_pxx'] = pxx[0]\n",
    "    results_dict['measured_pyy'] = pyy[0]\n",
    "    results_dict['measured_pzz'] = pzz[0]\n",
    "    results_dict['measured_pxy'] = pxy[0]\n",
    "    results_dict['measured_pxz'] = pxz[0]\n",
    "    results_dict['measured_pyz'] = pyz[0]\n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = relax_box(lammps_command, system, potential,\n",
    "                           mpi_command = mpi_command,\n",
    "                           p_xx = pressure_xx, \n",
    "                           p_yy = pressure_yy, \n",
    "                           p_zz = pressure_zz,\n",
    "                           strainrange = strainrange,\n",
    "                           tol=convergence_tol,\n",
    "                           diverge_scale=divergence_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'lx', 'ly', 'lz', 'xy', 'xz', 'yz', 'E_coh', 'measured_pxx', 'measured_pyy', 'measured_pzz', 'measured_pxy', 'measured_pxz', 'measured_pyz'])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __length_unit__ is the unit of length to display values in.\n",
    "- __energy_unit__ is the unit of energy to display values in.\n",
    "- __pressure_unit__ is the unit of pressure to display values in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_unit = 'angstrom'\n",
    "energy_unit = 'eV'\n",
    "pressure_unit = 'GPa'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2. Print Ecoh and lattice constants of relaxed ucell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ecoh = -4.449999998349075 eV\n",
      "a = 3.519999437540435 angstrom\n",
      "b = 3.5199994375404366 angstrom\n",
      "c = 3.519999437540417 angstrom\n",
      "alpha = 90.0\n",
      "beta =  90.0\n",
      "gamma = 90.0\n"
     ]
    }
   ],
   "source": [
    "print('Ecoh =', uc.get_in_units(results_dict['E_coh'], energy_unit), energy_unit)\n",
    "\n",
    "box = am.Box(lx=results_dict['lx'], ly=results_dict['ly'], lz=results_dict['lz'],\n",
    "             xy=results_dict['xy'], xz=results_dict['xz'], yz=results_dict['yz'])\n",
    "\n",
    "print('a =', uc.get_in_units(box.a / sizemults[0], length_unit), length_unit)\n",
    "print('b =', uc.get_in_units(box.b / sizemults[1], length_unit), length_unit) \n",
    "print('c =', uc.get_in_units(box.c / sizemults[2], length_unit), length_unit) \n",
    "print('alpha =', box.alpha)\n",
    "print('beta = ', box.beta)\n",
    "print('gamma =', box.gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3. Check final system pressures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pxx = 3.0498088437975e-12 GPa\n",
      "Pyy = 3.0488307763301005e-12 GPa\n",
      "Pzz = 2.1086605910971996e-12 GPa\n",
      "Pyz = -5.5323121222754e-16 GPa\n",
      "Pxz = -3.6252693770542003e-16 GPa\n",
      "Pxy = -1.3217127937177e-16 GPa\n"
     ]
    }
   ],
   "source": [
    "print('Pxx =', uc.get_in_units(results_dict['measured_pxx'], pressure_unit), pressure_unit)\n",
    "print('Pyy =', uc.get_in_units(results_dict['measured_pyy'], pressure_unit), pressure_unit)\n",
    "print('Pzz =', uc.get_in_units(results_dict['measured_pzz'], pressure_unit), pressure_unit)\n",
    "print('Pyz =', uc.get_in_units(results_dict['measured_pyz'], pressure_unit), pressure_unit)\n",
    "print('Pxz =', uc.get_in_units(results_dict['measured_pxz'], pressure_unit), pressure_unit)\n",
    "print('Pxy =', uc.get_in_units(results_dict['measured_pxy'], pressure_unit), pressure_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.4. Show relaxed atomic configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [10.560,  0.000,  0.000]\n",
      "bvect =  [ 0.000, 10.560,  0.000]\n",
      "cvect =  [ 0.000,  0.000, 10.560]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 108\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos', 'atom_id']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n",
      "      4 |       1 |   3.520 |   0.000 |   0.000\n",
      "      5 |       1 |   3.520 |   1.760 |   1.760\n",
      "      6 |       1 |   5.280 |   0.000 |   1.760\n",
      "      7 |       1 |   5.280 |   1.760 |   0.000\n",
      "      8 |       1 |   7.040 |   0.000 |   0.000\n",
      "      9 |       1 |   7.040 |   1.760 |   1.760\n",
      "     10 |       1 |   8.800 |   0.000 |   1.760\n",
      "     11 |       1 |   8.800 |   1.760 |   0.000\n",
      "     12 |       1 |   0.000 |   3.520 |   0.000\n",
      "     13 |       1 |   0.000 |   5.280 |   1.760\n",
      "     14 |       1 |   1.760 |   3.520 |   1.760\n",
      "     15 |       1 |   1.760 |   5.280 |   0.000\n",
      "     16 |       1 |   3.520 |   3.520 |   0.000\n",
      "     17 |       1 |   3.520 |   5.280 |   1.760\n",
      "     18 |       1 |   5.280 |   3.520 |   1.760\n",
      "     19 |       1 |   5.280 |   5.280 |   0.000\n",
      "     20 |       1 |   7.040 |   3.520 |   0.000\n",
      "     21 |       1 |   7.040 |   5.280 |   1.760\n",
      "     22 |       1 |   8.800 |   3.520 |   1.760\n",
      "     23 |       1 |   8.800 |   5.280 |   0.000\n",
      "     24 |       1 |   0.000 |   7.040 |   0.000\n",
      "     25 |       1 |   0.000 |   8.800 |   1.760\n",
      "     26 |       1 |   1.760 |   7.040 |   1.760\n",
      "     27 |       1 |   1.760 |   8.800 |   0.000\n",
      "     28 |       1 |   3.520 |   7.040 |   0.000\n",
      "     29 |       1 |   3.520 |   8.800 |   1.760\n",
      "     30 |       1 |   5.280 |   7.040 |   1.760\n",
      "     31 |       1 |   5.280 |   8.800 |   0.000\n",
      "     32 |       1 |   7.040 |   7.040 |   0.000\n",
      "     33 |       1 |   7.040 |   8.800 |   1.760\n",
      "     34 |       1 |   8.800 |   7.040 |   1.760\n",
      "     35 |       1 |   8.800 |   8.800 |   0.000\n",
      "     36 |       1 |   0.000 |   0.000 |   3.520\n",
      "     37 |       1 |   0.000 |   1.760 |   5.280\n",
      "     38 |       1 |   1.760 |   0.000 |   5.280\n",
      "     39 |       1 |   1.760 |   1.760 |   3.520\n",
      "     40 |       1 |   3.520 |   0.000 |   3.520\n",
      "     41 |       1 |   3.520 |   1.760 |   5.280\n",
      "     42 |       1 |   5.280 |   0.000 |   5.280\n",
      "     43 |       1 |   5.280 |   1.760 |   3.520\n",
      "     44 |       1 |   7.040 |   0.000 |   3.520\n",
      "     45 |       1 |   7.040 |   1.760 |   5.280\n",
      "     46 |       1 |   8.800 |   0.000 |   5.280\n",
      "     47 |       1 |   8.800 |   1.760 |   3.520\n",
      "     48 |       1 |   0.000 |   3.520 |   3.520\n",
      "     49 |       1 |   0.000 |   5.280 |   5.280\n",
      "     50 |       1 |   1.760 |   3.520 |   5.280\n",
      "     51 |       1 |   1.760 |   5.280 |   3.520\n",
      "     52 |       1 |   3.520 |   3.520 |   3.520\n",
      "     53 |       1 |   3.520 |   5.280 |   5.280\n",
      "     54 |       1 |   5.280 |   3.520 |   5.280\n",
      "     55 |       1 |   5.280 |   5.280 |   3.520\n",
      "     56 |       1 |   7.040 |   3.520 |   3.520\n",
      "     57 |       1 |   7.040 |   5.280 |   5.280\n",
      "     58 |       1 |   8.800 |   3.520 |   5.280\n",
      "     59 |       1 |   8.800 |   5.280 |   3.520\n",
      "     60 |       1 |   0.000 |   7.040 |   3.520\n",
      "     61 |       1 |   0.000 |   8.800 |   5.280\n",
      "     62 |       1 |   1.760 |   7.040 |   5.280\n",
      "     63 |       1 |   1.760 |   8.800 |   3.520\n",
      "     64 |       1 |   3.520 |   7.040 |   3.520\n",
      "     65 |       1 |   3.520 |   8.800 |   5.280\n",
      "     66 |       1 |   5.280 |   7.040 |   5.280\n",
      "     67 |       1 |   5.280 |   8.800 |   3.520\n",
      "     68 |       1 |   7.040 |   7.040 |   3.520\n",
      "     69 |       1 |   7.040 |   8.800 |   5.280\n",
      "     70 |       1 |   8.800 |   7.040 |   5.280\n",
      "     71 |       1 |   8.800 |   8.800 |   3.520\n",
      "     72 |       1 |   0.000 |   0.000 |   7.040\n",
      "     73 |       1 |   0.000 |   1.760 |   8.800\n",
      "     74 |       1 |   1.760 |   0.000 |   8.800\n",
      "     75 |       1 |   1.760 |   1.760 |   7.040\n",
      "     76 |       1 |   3.520 |   0.000 |   7.040\n",
      "     77 |       1 |   3.520 |   1.760 |   8.800\n",
      "     78 |       1 |   5.280 |   0.000 |   8.800\n",
      "     79 |       1 |   5.280 |   1.760 |   7.040\n",
      "     80 |       1 |   7.040 |   0.000 |   7.040\n",
      "     81 |       1 |   7.040 |   1.760 |   8.800\n",
      "     82 |       1 |   8.800 |   0.000 |   8.800\n",
      "     83 |       1 |   8.800 |   1.760 |   7.040\n",
      "     84 |       1 |   0.000 |   3.520 |   7.040\n",
      "     85 |       1 |   0.000 |   5.280 |   8.800\n",
      "     86 |       1 |   1.760 |   3.520 |   8.800\n",
      "     87 |       1 |   1.760 |   5.280 |   7.040\n",
      "     88 |       1 |   3.520 |   3.520 |   7.040\n",
      "     89 |       1 |   3.520 |   5.280 |   8.800\n",
      "     90 |       1 |   5.280 |   3.520 |   8.800\n",
      "     91 |       1 |   5.280 |   5.280 |   7.040\n",
      "     92 |       1 |   7.040 |   3.520 |   7.040\n",
      "     93 |       1 |   7.040 |   5.280 |   8.800\n",
      "     94 |       1 |   8.800 |   3.520 |   8.800\n",
      "     95 |       1 |   8.800 |   5.280 |   7.040\n",
      "     96 |       1 |   0.000 |   7.040 |   7.040\n",
      "     97 |       1 |   0.000 |   8.800 |   8.800\n",
      "     98 |       1 |   1.760 |   7.040 |   8.800\n",
      "     99 |       1 |   1.760 |   8.800 |   7.040\n",
      "    100 |       1 |   3.520 |   7.040 |   7.040\n",
      "    101 |       1 |   3.520 |   8.800 |   8.800\n",
      "    102 |       1 |   5.280 |   7.040 |   8.800\n",
      "    103 |       1 |   5.280 |   8.800 |   7.040\n",
      "    104 |       1 |   7.040 |   7.040 |   7.040\n",
      "    105 |       1 |   7.040 |   8.800 |   8.800\n",
      "    106 |       1 |   8.800 |   7.040 |   8.800\n",
      "    107 |       1 |   8.800 |   8.800 |   7.040\n"
     ]
    }
   ],
   "source": [
    "finalsystem = am.load('atom_dump', results_dict['dumpfile_final'],\n",
    "                      symbols=results_dict['symbols_final'])\n",
    "print(finalsystem)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
